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Abstract. In this note we describe a procedure of calculating the number all regular tetrahedra 
that have coordinates in the set {0, 1, n}. We develop a few results that may help in finding good 
estimates for this sequence which is twice A103158 in the Online Encyclopedia of Integer Sequences 



1. INTRODUCTION 

The story of regular tetrahedra having vertices of integer coordinates starts with the parametriza- 
tion of some equilateral triangles in Z 3 that begun in There was an additional hypothesis that 
did not cover all the generality in the result obtained in [9] but it was removed successfully in [2]. 
A few other related results appeared in [10] and [11]. In this note we are interested in the following 
problem. 

How many regular tetrahedra, T(n), can be found if the coordinates of its vertices must be in the 
set {0, 1, ...,n}? We observe that ^4103158 = ^T(n), see [13J. This sequence starts as in the next 
tables: 



n 


1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


11 


A103158 


1 


9 


36 


104 


257 


549 


1058 


1896 


3199 


5154 


7926 



X 



n 


12 


13 


14 


15 


16 


17 


18 


A103158 


11768 


16967 


23859 


32846 


44378 


58977 


77215 



Using our method which is going to be described later we extended this sequence for all n < 100, 
and one can go far enough with this if time allows and powerful computer is used. 

The rest of the terms are included at the end of the paper. Our approach begins with looking first 
at the faces of a regular tetrahedron, which must be equilateral triangles. It turns out that every 
equilateral triangle in Z 3 after a translation by a vector with integer coordinates can be assumed 
to have the origin as one of its vertices. Then one can show that the other triangle's vertices are 
contained in a lattice of points of the form 
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(1) V aAc :={{a,(3,~f) eZ 3 \ aa + b{3 + = 0, a 2 + b 2 + c 2 = 3d 2 , a,b,c,deZ}. 




Figure 1: The lattice V a .b. c 

In general, the vertices of the equilateral triangles that dwell in V a {, c form a strict sub-lattice of 
^a,6,c which is generated by only two vectors, £ and if (see Figure [I]). These two vectors are 
described by the Theorem 11.11 proved in [2]. 

Theorem 1.1. Let a, b, c, d be odd integers such that a 2 + b 2 + c 2 = 3d 2 and gcd(a, b, c) = 1. 
Then for every m,n £ 7L (not both zero) the triangle OPQ, determined by 

(2) OP = mC,-nrf, OQ = n ( - (n - m)lf , with ( = (Ci, Ci , C2), ~V = (Vi, V2, m), 



(3) 



rac + dbs 



das — bcr 



X3 = r, 



db(s — 3r) + ac(r + s) 
da(s — 3r) — bc(r + s) 



2q 



r + s 



where q = a 2 + b 2 and (r, s) is a suitable solution of 2q = s 2 + 3r 2 that makes all the numbers in 
([3j) integers, forms an equilateral triangle in Z 3 contained in the lattice ([1]) and having sides-lengths 



equal to dy / 2(m 2 — mn + n 2 ). 
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Conversely, there exists a choice of the integers r and s such that given an arbitrary equilateral 
triangle in M 3 whose vertices, one at the origin and the other two in the lattice ([1]), then there also 
exist integers m and n such that the two vertices not at the origin are given by © and © . 

The Diophantine equation 



(4) 



a 2 + b 2 + c 2 



3d 2 



has non-trivial solutions for every d odd. As a curiosity, for d = 2009 one obtains 294 solutions 
satisfying also < a < b < c and gcd(a, b, c) = 1. We will refer to such a solution of (J2J) as 
a positive ordered primitive solution. For d = 2008, all these solutions satisfy even a stronger 
condition: a < b < c. Determining the exact number of solutions for Q is certainly important 
if one wishes to find the number (or just an estimate) of equilateral triangles or the number of 
tetrahedra with vertices in {0, 1, 2, n} 3 . The number of solutions for Q, coincidentally, taken 
into account all permutations and changes of signs is given in a 1999 paper of Hirschhorn and Seller 



(5) 



n 

p = 1 or 7(mod 12) 
pP\\d 



n 

5 or ll(mod 12) 
q a \\d 



q a + 2- 



{1 % J~ 3 1 d 
„ 7 _ 1 , Even more important for our purpose is the calculation of the 

2 if | \ d. 

number of primitive representations of d as in (gcd(a, b, c) = 1) in terms of d which appeared 
in a very recent paper of Cooper and Hirschhorn [3j. One may check easily that the following is a 
corollary of Theorem 2 in [3] . 

Theorem 1.2. [Cooper-Hirschhorn] Given an odd number d, the number of primitive so- 
lutions of dU) taking into account all changing of sings and permutations, is equal to 



(6) 



A(d) := 8d ]J ( 1 

p\d,p prime \ 



P 



where (-^) it the Legendre symbol. 



We remind the reader that, if p is prime then 
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(7) 



(^) 



ifp = 3 

1 if p = 1 or 7 (mod 12) 



-1 if p = 5 or 11 (mod 12) 

We observe that the same type of prime partition are used into different calculations in both 
formulae ([5]) and ([6]). We have mentioned that the number of positive ordered primitive rep- 
resentations for d = 2009 was 294. This is exactly the number given by ([6]) modulo the num- 
ber of permutations and changes of signs: indeed, 2009 = (41)(7 2 ), {-f-) = 1, (-gj) = —1 and 
A ^g 09 ^ = 8 ^g 7 - (1 — ^)(1 + ji) = 7(42) = 294. This happens because there are no repeating 
values for a, b and c in any of the positive ordered primitive solutions of We will see later how 
the correct number or positive ordered primitive representations can be obtained in general. 



For k G N, we let Q. := {(m,n) 6ZxZ: m 2 



■mn + n 



k 2 }. In |9j we showed that every regular 



tetrahedron with integer coordinates must have side lengths of the form A\/2, A G N, and in [TT] 
we have found the following characterization of the regular tetrahedrons with integer coordinates. 

Theorem 1.3. Every tetrahedron whose side lengths are Av2, A G N, which has a vertex at 
the origin, can be obtained by taking as one of its faces an equilateral triangle having the origin as 
a vertex and the other two vertices given by ([2]) and ([3]) with a, b, c and d odd integers satisfying 

with d a divisor of A, and then completing it with the fourth vertex R with coordinates 



(8) 



/ (2Ci - Vi) m 

-(Ci + m) n 

±2ak 



(2( 2 - m) m 
-{C,2 + m)n 
±2bk 



(2C3 - m) m \ 
±2ck 



for some (m,n) G O(fc), k :- 



A 

7r 



V / 

Conversely, if we let a, b, c and d be a primitive solution o/dH), let k G N and (m,n) G Q(k), then 
the coordinates of the point R in ([ED, which completes the equilateral triangle OPQ given as in (|2|) 
and ([3]), are 

(a) all integers, if k = (mod 3) regardless of the choice of signs or 

(b) integers, precisely for only one choice of the signs if k (mod 3). 

The following graph (Figure 2) is constructed on the positive ordered primitive solutions of ((4]), 
with edges defined by: 

two vertices, say [(oi, &i, ci), d{\ and [(02, C2), dz], are connected, if and only if 



COUNTING ALL REGULAR TETRAHEDRA IN {0,1,..., n} 3 



5 



0) 



a\a 2 ± a2&2 ± cic 2 ± dic?2 = 



for some choice of the signs and permutation (a^b^, c' 2 ) of (02, b 2 , c 2 ). 



Equation ([9]) insures basically that the planes "P aij b ljCl and V, 



a 2>^2' c 2 



associated to two faces make 



a dihedral angle of arccos(l/3) ~ 70.52878°. In fact, this equality characterizes the existence of 

a regular tetrahedron having integer coordinates with one of its faces in the plane 7 3 aii 6 liCl and 

another contained in the plane V a > t b' ,c' 2 - For instance, [(1, 1, 5), 3] is connected to [(1, 5, 11), 7] since 

1(11) + (1)5 + 5(1) — 3(7) = 0. An example of a regular tetrahedron which has a face in V-5-i t i 

and one face in V-i,-&,u is given by the vertices: [19, 23, 0], [0, 12, 20], [27, 0, 17], and [24, 27, 29]. 

[(1,11,31),19] 

[(13,17,25),1 9] [(1,5,29),17] 

R7.17.23).17l \ 

[(11 ,11 ,25^1 7] 
[(13,1 3,23), 17] 




[(1.1.1). 



,11,19),13] 
j(7,13,17),13] 

[(1,1,19),11] 

[(5,7,17),11] 
[(5,'13,13),11] 



[(1.1.5),3] 



[(5,7,1 3),9] 



Figure 2: The graph TZT , d < 19. 

A few questions related to this graph appear naturally at this point. Is it connected? Is there 
a different characterization of the existence of an edge between two vertices in terms of only d\ 
and d 2 ? We do not have an answer to the second question. This graph seems to have a fractal 
structure. 
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Each edge in this graph, determined by [(a±, bi, c\), d±] and [(02, 62, C2), tfe], gives rise to a minimal 
tetrahedra (the side lengths are at most max{di, cfojv^) which is determined up to the set of 
isometric transformations that are generated by the symmetries of the cube C(m) where m is the 
size of the smallest "cube" {0, 1, • • • , m} 3 containing the tetrahedron or a translation of it. 

2. Some preliminaries 

We would like to have a good estimate of the primitive solutions of (jlj) which satisfy in addition 
< a < b < c. Let us observe that we cannot have a = b = c unless d = 1. So, the counting in @ 
via (0) would give what we want if we can count the number of positive primitive solutions of the 
following equation in terms of d: 



(10) 2a 2 + c 2 = 3d 2 . 

A similar description to the Pythagorean triples, which gives the nature of the solutions of (|10j) . is 
stated next. 

Theorem 2.1. For every positive integers I and k such that, gcd(k,l) = 1 and k is odd, then 
a, c and d given by 



(11) d = 2l 2 + k 2 and < 



\2l 2 + 2kl - k% c=\k 2 + 4kl - 2l\ if k^l (mod 3) 



a = \2l 2 - 2kl - k 2 \, c=\k 2 - 4kl - 2l 2 \, if k^-l (mod 3) 
constitute a positive primitive solution for (|10p . 

Conversely, with the exception of the trivial solution a = c = d = 1, every positive primitive 
solution for (|10p appears in the way described above for some I and k. 



Proof. First, one can check that (llip satisfy (llOh for every I and k. As a result it follows 
that a, c and d are positive integers. Let p be a prime dividing a, c and d. Then p must divide 
±a — d = 2k(±l — k) and so p is equal to 2, p divides k or it divides ±Z — k. Yip = 2 then, p must 
divide k but this contradicts the assumption that k is odd. 

In case p is not equal to 2 and it divides k, we see p must divide I 2 = (d — k 2 )/2. Since we 
assumed gcd(l, k) = 1 it remains that p must divide ±Z — k. By our assumptions on k and I, p 
cannot be equal to 3. Then p divides ±a + (±Z — k) 2 = 3Z 2 . Because p 7^ 3 then p must divide Z 2 
and so p should divide Z and then k. This contradiction shows that a, c and d cannot have prime 
common factors. So, we have a primitive solution in (jlip . 

For the converse, let us assume that a, c and d is a positive primitive solution of (|10p . which is 
different of the trivial one. We denote by u = % and v = § . Then the point of rational coordinates 
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(u, v) (different of (1, 1))) is on the ellipse jj^ + ^- = 1 (Figure 3) in the first quadrant. This ellipse 
contains the following four points with integer coordinates: (1, 1), ( — 1, 1), ( — 1, —1) and (1, 1). This 
gives the lines y + 1 = t\{x + 1), y + 1 = t 2 {x — 1), y — 1 = tz(x + 1), and y — 1 = i4(x — 1), 
passing through (ii,t>) and one of the points mentioned above. Hence, the slopes t%, t 2 , £3, and £4, 
are rational numbers. This gives expressions for the point (u, v) in terms of ti {i = 1, ...,4). Let us 
assume that i« = t 1 with ki, U 6 Z, written in the reduced form. Then we must have 

_ \2± 2tj -tj\ _ \2lf ± 2kjlj - kf\ _ \tj ± 4tj - 2| _ |fc, 2 ± 4fc^ - 2lf\ 
U ~ 2 + t} " 2fTfcf ' V ~ 2 + t} " 2fTfcf ' 




Figure 3: The ellipse + \ = 1 

and so, these equalities give 

, 10 . a \2%±2kili-k?\ ,c \k$±4kih-2l?\ . 1 , 

(12) d = 2i> + k? ' and 3 = 2fTfef ' * = h -' 4 - 

We claim that the function t{ — > 2lf + kf (i = 1, 4) is injective. If for some 2Z| + kf = 21? + fcj 
(* 7^ i)> that would imply that the corresponding numerators in (|12p are equal. This gives enough 
infirmation to conclude a contradiction. There are („) = 6 possibilities here but we are going to 
include the details only in the case i = 1 and j = 2. The rest of the cases can be done in a similar 
fashion. For this situation we have, 2l\ + 2k\l\ — k\ = k\ + 2^2 — 2l\ and k\ + 4k\li — 2l\ = 
k\ — 4k 2 l 2 — 2l\ . The first equality implies 

2k x l\ = k\ + kl + 2k 2 l 2 - 2l\ - 2l\ = 2k\ + 2k 2 l 2 - 
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which substituted into the second equality gives 



6kf = 2k\ - 8k 2 l 2 + 81% O 3fc? = (k 2 - 2l 2 ) 2 . 
Because \/3 is irrational, the last equality is impossible for k\, k 2 , l 2 integers and k\ nonzero. 
For the other cases one will get a contradiction based on the facts that y| and \[2 are irrational 
numbers. 

A similar argument to the one in the first part of the proof shows that the fractions in the right- 
hand side of the equalities of (|12p can be simplified only by a factor of 2, 3 or 6. Having four distinct 
possibilities in (|12p for the denominators, exactly one of the fractions (simultaneously in the first 
and second equalities) must be in reduced form. This one will give the wanted representation. ■ 

Similar to Fermat's theorem about the representation of primes as a sum of two squares and the 
number of such representations one can show the next result. 



(u) r 2 (d) = { 



Theorem 2.2. (Fermat [3]) An odd prime p can be written as 2x 2 + y 2 with x,y G Z if and 
only ifp = l or 3 (mod 8). If d = 2 k ]]p e * i ]\qf is the prime factorization of d, with qj primes 
as before and pi the rest of them, then the number of representations d = 2x 2 + y 2 with x, y G Z is 
either zero if not all o.{ are even and otherwise given by 

as) L^n^ +i )j- 

The number of positive primitive representations d = 2x 2 + y 2 for d odd, i.e. x,y G N and 
gcd(x,y) = 1, is equal to 

if d is divisible by a prime factor of the form 8s+5 or 8s+7, s > 0, 
where k is the number of distinct prime factors of d 

2 fc-l 

of d of the form 8s+l, or 8s+3 (s > 0). 

Putting the two results together (Theorem 12.21 and Theorem II. 2p we obtain the following propo- 
sition: 

Proposition 2.3. For every d odd, the number of representations of which satisfy < a < 
b < c and gcd(a, b, c) = 1 is equal to 

A(d) + 24r 2 (3d 2 ) 

(15) vre(d) = — . 

A regular tetrahedron whose vertices are integers is said to be irreducible if it cannot be obtained 
by an integer dilation and a translation from a smaller one also with integer coordinates. An 
important question at this point about irreducible tetrahedra is included next. 
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Does every irreducible tetrahedron with integer coordinates have a face with a normal 
vector (a,b,c) satisfying a 2 + b 2 + c 2 = 3d 2 , such that d gives the side lengths I of 
the tetrahedron by the formula £ = dy/2? In other words, is there a face for which 
k = 1 in the Theorem \1.3\ ? 

Unfortunately the answer to this question is no. The following points together with the origin, 
[-6677, -2672, 1445], [-5940, 4143, -1167], [-3837, 2595, 5688] form a regular tetrahedron of side- 
lengths equal to 5187\/2 and the highest d for the faces is 1729. We observe that 3, 7, 13 and 19 
are the first three distinct primes numbers of the form u 2 + 3v 2 , «,!)6Z. 



The program is written in Maple code and it is based on the Theorem 11.31 The main idea is 
to create a list of irreducible regular tetrahedra that can be used to generate all the others in 
{0, 1, 2, n} 3 by certain transformations generating a partition for the set of all the tetrahedra. 
Each such irreducible tetrahedron is constructed out of the equation of one face using Theorem ll.31 
One important problem that appears here is to make sure this list contains distinct elements, 
elements which may appear theoretically in this list from four different constructions, one for each 
face. It turns out that there is a simple way of making sure that this doesn't happen. 

Let us observe that if gcd(m, n) = d > 1 then all the coordinates of the vertices of the initial 
face are multiple of d and by the formula (|8|), so are the coordinates of the fourth point of the 
tetrahedron. This is because the numbers k in (JHJ) satisfy a Diophantine equation of the form 
k 2 = m 2 — mn + n 2 . We go one step further, if k 2 is of the form 



then one can see that for k even, the gcd(m, n) > 2. Also, if k is odd then both m and n must be 
odd and so we have a representation of k 2 as u 2 + 3v 2 , u,v G Z. If k is divisible by 3 then it is easy 
to see that 3 divides u an v, and this attracts gcd(m, n) > 3. Hence we are going to look only for 
those odd values k < n, which are not multiples of 3, in the Theorem 11.31 This means that only 
one choice of signs in (JSj) is useful. A similar fact to Theorem 12.21 takes place. 

Theorem 3.1. (Fermat [4J) A prime p can be written as x 2 + 3y 2 with x,y £ 7L if and only if 
p = 3 or p = 1 (mod 3). If d = Ylp^ Y\ Qj j with qj primes as before and p. L the rest of them, then 
the number of representations d = x 2 + 3y 2 with x,y G 7L is either zero if not all ai are even and 
otherwise given by 



3. The Code 



(16) 
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(17) 



The number of positive primitive representations d 
gcd(x,y) = 1, is equal to 



x 2 + 3y 2 for d odd, i.e. x,y G N and 



if d is even or divisible by a prime factor of the form 3s+2, s > 



(is) r 3 (d) 



< 



where k is the number of distinct prime factors of d such as 3 or 



2 



< 



[ of the form 3s+l (s > 2). 

As a result of these facts we first calculate all k < n such that k 2 = m 2 — ran + n 2 has a solution 
with gcd(m,n) = 1. 

kvalues:=pioc(n) 



local i,j,k,L,a,p,q,r,m,mm; 
L:=;mm:=floor((n+l) /2); 

for i from 2 to mm do a:=ifactors(2i-l); k:=nops(a[2]);r:=0; 
for j from 1 to k do 

m:=a[2] [j] [1]; p:=m mod 3; 

if m=3 then r:=l fi; 

if r=0 and p=2 then r:=l fi; 



The result of this procedure for n = 100 is [1, 7, 13, 19, 31, 37, 43, 49, 61, 67, 73, 79, 91, 97]. These 
are all the natural numbers less than 100 which are primes of the form 3s + 1 or products of such 
primes. We noticed that it takes only a fraction of a second to execute this procedure if we limit n 
to be less than 10000, although it may time consuming for big numbers. So an alternative solution 
to this procedure may use a similar result to that in Theorem 12.11 to describe all the solutions 
of (|16p . For each k found by the previous procedure, there are usually at least eighteen solutions 
of (|16p if signs and order are counted, but if we impose the conditions gcd(m,n) = 1, < m, n 
and 2m < n, we slice these solutions by a factor of 18. Such a solution is going to be referred to 



od; 

if r=0 then L:=L union {2i — l};fi; 
od; 

L:=L union {1}; L:=convert(L,list); 
end: 
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as a primitive solution and these primitive solution of (fl~6j) can be calculated with the following 
procedure. 

listofmn:=proc(k) 
local a,b,i,x,m,n,nx,L; 

x:=[isolve(/c 2 = m? — ran + n 2 )]; nx:=nops(x); L:={ }; 
for i from 1 to nx do 
if lhs(x[i][l])=m then a:=rhs(x[i] [1]); b:=rhs(x[i] [2]); 

else b:=rhs(x[i][l]); a:=rhs(x[i] [2]); fi; 
if gcd(a,b)=l and a > and b > and 2a < b then L:=L union [a,b];fi; 
od;L; 
end: 

For example, if k = 91 we get two primitive solutions, m = 1991, n = 9095, and m = 3401, 
n = 9440. It turns out that it is enough to know just the primitive solutions of (|16p in order to find 
all integer solutions. For each k < n, output of the procedure kvalues, we need to find the values of 
d as in Theorem ! 1.3} which are only restricted to two conditions: d must be an odd positive integer 
and dk < n. The last restriction follows from the fact that the sides of the regular tetrahedron 
constructed from say d and k, with a primitive solution (m, n) of (|16p . must be equal to dk\[2. 
One can see that the biggest regular tetrahedron inscribed in the cube [0, n] 3 has sides equal to 
(see Proposition 2.1 in [10]). The next procedure is then very simple. 



D C 




Figure 4: Largest tetrahedrons in a cube 

determined:=proc(n) 
local i,x,m,L,j; 
x: =kvalues (n) ;m: =nops (x) ; 
for i from 1 to m do 
j:=l;L[i]:={}; 
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while (2j-l)x[i]<=n do 

L[i]:=L[i] union {2j - l};j:=j+l; od; 
od; [seq(L[i],i=l..m)]; 
end: 

Next, we need to find all primitive solutions (a, b, c) of the equation a 2 + b 2 + c 2 = 3d 2 . The number 
of such solution is given by ([To]) . 

abcsol: =proc(d) 

local i,j,k,m,u,x,y,sol,cd; 

sol:={}; 

for i from 1 to d do 
u:=[isolve(3d 2 — i 2 = x 2 + y 2 )];k:=nops(u); 
for j from 1 to k do 

if rhs(u[j][l])>=i and rhs(u[j] [2])>=i then 
cd: =gcd(gcd(i,rhs(u[j] [1] ) ) ,rhs(u [j] [2] ) ) ; 

if cd=l then sol:=sol union {sort([i, rhs(u[j][l]), rhs(u[j][2])])};R; 

fi; 
od; od; 

convert (sol, list); 
end: 

For example, if d = 2009 we get 294 solutions, as seen before, of which one of them is a = 1, 
b = 1159 and c = 3281. Next, we are going to use the construction of an equilateral triangle in 
the plane of equation ax + by + cz = with a 2 + b 2 + c 2 = 3d 2 using the formulae ([3|) with m and 
n given by the procedure listofmn. The fourth point is then calculated using the formula ([8]). We 
are using only two of the possible values of m and n {m 1 = m, n' = n and m' = n, n' = n — m, 
in order to obtain two equilateral triangles that share a side, OQ) since all other tetrahedra as in 
the figure below, can be obtained from these two by a simple translation, and as a result they will 
translate into the same minimal tetrahedron inside the first quadrant (see |llj). 
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Figure 5: All six tetrahedrons with one face in V a ,b,c 
/mdpar:=proc(a,b,c,mm,nn) 

local ij,sol,mx,nx,r,s,my,ny,q,d,u,v,w,x,y,z,mu,nu,mv,nv,ef,ns,mz,nz, 

mw,nw,oml,om2,l,uu,t,Rl,R2,fc,k; 

q := a 2 + b 2 ; k := sqrt(mm 2 — mm * nn + nil 2 ); 

sol:=convert(isoive(2*g = x 2 + 3*y 2 ),list); ns:=nops(sol); d:=sqrt((a 2 + 6 2 + c 2 )/3); 
ef:=0; 

for i from 1 to ns do 

if ef=0 then r:=rhs(sol[i][l]); s:=rhs(sol[i][2]); 
uu: = (s 2 + 3 * r 2 — 2 * q) 2 ; if uu > then t:=s;s:=r;r:=t; fi; 
mx:=-(d *b* (3*r + s) + a* c* (r — s))/(2* g);nx:=-(r * a * c + d * b * s) / q; 
my: = (d *a*(3*r + s) — b* c* (r — s))/(2 * g);ny:=-(r * b* c — d* a* s) jq\ 
mz:=(r-s)/2;nz:=r; mu:=nx;mv:=ny;mw:=nz;nu:=nx-mx;nv:=riy-my; nw:=nz-mz; 
if mx=floor(mx) and nx=floor(nx) and my=floor(my) and ny=floor(ny) then 
u:=mu*m-nu*n;v:=mv*m-nv*n;w:=mw*m-nw*n; 
x:=mx*m-nx*n;y:=my*m-ny*n;z:=mz*m-nz*n; 
Rl:=[(x+u-2*a*k)/3,(v+y-2*b*k)/3,(z+w-2*c*k)/3]; 
R2:= [(x+u+2*a*k) /3,(v+y+2*b*k) /3, (z+w+2*c*k) /3] ; 
fc:=subs(m=mm,n=nn,Rl[l]); if fc=floor(fc) then 
oml:=subs(m=mm,n=nn,[[u,v,w],[x,y,z],Rl]); else 
oml:=subs(m=mm,n=nn,[[u,v,w],[x,y,z],R2]);fi; 
fc:=subs(m=nn,n=nn-mm,Rl[l]); if fc=floor(fc) then 
om2:=subs(m=nn,n=nn-mm,[[u,v,w],[x,y,z],Rl]); else 
om2:=subs(m=nn,n=nn-mm,[[u,v,w],[x,y,z],R2]);fi; ef:=l; fi;fi; od; 
oml,om2; end: 
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Before we translate these two tetrahedrons we need a small sub-routine for subtraction of two 
vectors. 

subtrv:=proc(XJ ,V) local W; 
W[l] :=U[1]-V[1] ;W[2] :=U[2]-V[2] ;W[3] :=U[3]-V[3] ; [W[l] ,W[2] ,W[3]] ; end: 
The next procedure translates a tetrahedron which comes as output of findpar into the positive 
octant of the space in such a way that for each component at least one of the tetrahedron's vertex has 
a zero coordinate on that component. Let us observe that this operation is invariant to translations 
of the tetrahedron; this justifies the choice of looking at only two tetrahedrons out of six in the 
procedure findpar. 

tmttopq:=\)Toc(T) local i,a,b,c,v,0,TR; 
a:=min(T[l][l],T[2][l],T[3][l],0); 
b:=min(T[l][2],T[2][2],T[3][2],0); 

c:=min(T[l][3],T[2][3],T[3][3],0); O:=[0,0,0]; v:=[a,b,c]; 
TR:={subtrv(0,v),subtrv(T[l],v),subtrv(T[2],v),subtrv(T[3],v)}; end: 

Next, we calculate the size of smallest cube C m := [0, to] 3 which contains the tetrahedron resulted 
from the tmttopq. 

msco/mt:=proc(Q) local a,b,c,T,m; T:=convert(Q,list); 
a:=max(T[l][l],T[2][l],T[3][l],T[4][l]); 
b:=max(T[l][2],T[2][2],T[3][2],T[4][2]); 

c:=max(T[l][3],T[2][3],T[3][3],T[4][3]); m:=max(a,b,c); end: 
The tetrahedron obtained as a result of tmttopq is then transformed within the cube found above 
through all the translations, rotations and symmetries of the cube. We denote this orbit of T, by 
0(T). 

orbitl:=proc(T) local i,k,Tl,a,b,c,x,T2,T3,T4,T5,T6,T7,T8,T9,T10,Tll,T12,T13,T14,T15,T16, 
T17,T18,T19,T20,T21,T22,T23,T24,S,Q,d,al,bl,cl; Q:=convert(T,list); 
d:=mscofmt(T); T1:=T; T2:={seq([Q[k][2],Q[k][3],Q[k][l]],k=1..4)}; 
T3:={seq([Q[k][l],Q[k][3],Q[k][2]],k=1..4)}; 
T4:={seq([Q[k][l],Q[k][2],d-Q[k][3]],k=1..4)} 
T5:={seq([Q[k][2],Q[k][3],d-Q[k][l]],k=1..4)} 
T6:={seq([Q[k][l],Q[k][3],d-Q[k][2]],k=1..4)} 
T7:={seq([Q[k][l],d-Q[k][2],Q[k][3]],k=1..4)} 
T8:={seq([Q[k][2],d-Q[k][3],Q[k][l]],k=1..4)} 
T9:={seq([Q[k][l],d-Q[k][3],Q[k][2]],k=1..4)} 
T10:={seq([d-Q[k][l],Q[k][2],Q[k][3]],k=1..4)}; 
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Til 
T12 
T13 
T14 
T15 
T16 
T17 
T18 
T19 
T20 
T21 
T22 
T23 
T24 



={seq( 
={seq( 
={seq( 
={seq( 
={seq( 
={seq( 
={seq( 
={seq( 
={seq( 
={seq( 
={seq( 
={seq( 
={seq( 
={seq( 



;d-Q[k][2],Q[k][3],Q[k][l]],k=1..4)}; 
;d-Q[k][l],Q[k][3],Q[k][2]],k=1..4)}; 



[Q[k][l],d-Q[k][2],d-Q[k][3]],k=1..4)} 
[Q[k][2],d-Q[k][3],d-Q[k][l]],k=1..4)} 
[Q[k][l],d-Q[k][3],d-Q[k][2]],k=1..4)} 
[d-Q[k][l],d-Q[k][2],Q[k][3]],k=1..4)} 
[d-Q[k][2],d-Q[k][3],Q[k][l]],k=1..4)} 
[d-Q[k][l],d-Q[k][3],Q[k][2]],k=1..4)} 
[d-Q[k][l],Q[k][2],d-Q[k][3]],k=1..4)} 
[d-Q[k][2],Q[k][3],d-Q[k][l]],k=1..4)} 
[d-Q[k][l],Q[k][3],d-Q[k][2]],k=1..4)} 
[d-Q[k][l],d-Q[k][2],d-Q[k][3]],k=1..4)} 
[d-Q[k][2],d-Q[k][3],d-Q[k][l]],k=1..4)| 
[d-Q[k][l],d-Q[k][3],d-Q[k][2]],k=1..4)} 
S:={T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11,T12,T13,T14,T15,T16,T17, 
T18,T19,T20,T21,T22,T23,T24};S; end: 
orbit:=proc(T) local S,Q,T1; Q:=convert(T,list); 
T1:={[Q[1][3],Q[1][2],Q[1][1]],[Q[2][3],Q[2][2],Q[2][1]], 
[Q[3][3],Q[3][2],Q[3][1]],[Q[4][3],Q[4][2],Q[4][1]]}; 
S:=orbitl(T) union orbitl(Tl); S; end: 

We recall from |10j a few variables that we are going to use in this calculation also. The theorem 
used there applies as well to this case because it is a pure set theoretic result. The meaning of 
those variables here is: 

(i) n -the dimension of the cube C n = [0, n] 3 , 

(ii) m -the maximum of all the coordinates in a tetrahedron T computed by tmttopq, 
(hi) a(T) -the cardinality of 0(T) within C m , 

(iv) (3(T) -the cardinality of 0(T) n [O(T) + ej, 

(v) 7 (T) -the cardinality of [0{T) + ei] n [0(T) + e 2 ]. 

Theorem 3.2. (Theorem 2.2 inflOjj The number f(T, n) of all tetrahedrons that can be obtained 
from T within a cube C n by translations, rotations, or symmetries, is given by 



(19) f(T,n) = (n + l-mfa-3{n + l 
for all n> m. 



m 



) 2 (n — m)/3 + 3(n + 1 — m)(n — m) 2 j, 



Hence, we need to calculate a, (3 and 7. 
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addvec: =proc(U,V) local W; 
W[1]:=U[1]+V[1];W[2]:=U[2]+V[2];W[3]:=U[3]+V[3];[W[1],W[2],W[3]]; 
end: 

addvect: =proc(T,v) 
local i,Q;Q:={ }; 
for i from 1 to 4 do 
Q:=Q union addvec (T[i],v); 
od;Q;end: 

transl:=pioc(T) 
local S,Q,i,j,k,a2,b2,c2,a,b,c,d; 
Q : = convert ( T , list ) ; 

a:=max(Q[l][l],Q[2][l],Q[3][l],Q[4][l]); 

b:=max(Q[l][2],Q[2][2],Q[3][2],Q[4][2]); 

c:=max(Q[l][3],Q[2][3],Q[3][3],Q[4][3]); 

d:=max(a,b,c); 

a2:=d-a;b2:=d-b;c2:=d-c; 

S:=orbit(T); 

for i from to a2 do 

for j from to b2 do 

for k from to c2 do 
S:=S union orbit (add vect(T,[i,j,k])); 
od; od; od; S; end: 

This last procedure gives the value of a. Then (3 is calculated by the following. 

inters along El: =proc(T) 
local S,m,i,Sl,S2; 

S2:=transl(T);S:=convert(S2,list);m:=nops(S);Sl:={ }; 
for i from 1 to m do S1:=S1 union addvect(S[i], [0,0,1]); 
od; S2 intersect SI; end: 

Then the variable 7 is given by the procedure: 

intersalongE2:=proc(T) local S,m,i,Sl,S2,S3,S4; 
S2:=transl(T);S:=convert(S2,list);m:=nops(S);Sl:={ }; 
for i from 1 to m do S1:=S1 union {addvect(S[i], [0,0,1])}; 
od; S3:={}; 
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for i from 1 to m do S3:=S3 union {addvect(S[i], [0,1,0])}; od; SI intersect S3; 
end: 

The function in Theorem 13.21 is then implemented by 

/:=(n,d, alpha, beta, gamma)^ (n — d + l) 3 * alpha — 3 * (n — d) * (— d + 1 + n) 2 * 
beta + 3 * gamma * (n — d + 1) * (n — d) 2 : 

notetraincn:=proc(T ,n) 
local d,x,y,z,w; 
d:=mscofmt(T); 

x:=nops(transl(T));y:=nops(intersalongEl(T));w:=nops(intersalongE2(T)); 
if n > d then 
z:=f(n,d,x,y,w);else z:=0; 
fi;z;end: 

Finally we put together a list of irreducible tetrahedrons whose orbits under the operations above 
are pairwise disjoint. The next four procedures are pretty simple and one can figure out what they 
do. They are used in the code of ExtendList. 

cfetance:=proc(A,B) local C; 

C:=subtrv(A,B); sqrt(C[l] 2 + C[2] 2 + C[3} 2 ); 
end: 

checkrt: =proc(T) local dl,d2,d3,d4,d5,d6,Q,D; 
dl:=distance(T[4],T[l])/sqrt(2); d2:=distance(T[4],T[2])/sqrt(2); 
d3:=distance(T[4],T[3])/sqrt(2); d4:=distance(T[l],T[2])/sqrt(2); 
d5:=distance(T[l],T[3])/sqrt(2); d6:=distance(T[2],T[3])/sqrt(2); 
min(dl,d2,d3,d4,d5,d6); end: 

crossproductt:=proc(\J ,V) local x,i,j,k,d,il,jl,kl,dl; 

i:=U[2]*V[3]-U[3]*V[2];j:=U[3]*V[l]-U[l]*V[3];k:=U[l]*V[2]-V[l]*U[2]; 

d:=gcd(i,j);dl:=gcd(d,k);il:=i/dl;jl:=j/dl;kl:=k/dl; sqrt((il 2 + jl 2 + £d 2 )/3); 
end: 

facesnew: =proc(T) local N1,N2,N3,N4,U,V,W; 
U:=subtrv(T[l],T[2]);V:=subtrv(T[l],T[3]);W:=subtrv(T[l],T[4]); 
Nl:=crossproductt(U,V);N2:=crossproductt(V,W);N3:=crossproductt(U,W); 
N4:=crossproductt(subtrv(U,V),subtrv(U,W)); max(Nl,N2,N3,N4); end: 

ExtendList:=proc(n,L,mm,im) 

local i, sol, nsol,alpha,beta,gammma,nel,mt,ttpl,ttp2, Or 1, Or, tnel,NL,intcard,cio, normals, length, k; 
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nel:=nops(L); 

k:=sqrt(mm 2 — mm *nn + nn 2 ); 
sol:=abcsol(n);nsol:=nops(sol); 
tnel:=nel;Or:={ }; 
for i from 1 to nsol do 

mt:=findpar(sol[i][l],sol[i][2],sol[i][3],mm,nn); ttpl:=tmttopq(mt[l]); ttp2:=tmttopq(mt[2]); 

normals: =facesnew (ttp 1 ) ;length: =checkrt (ftp 1 ) ; 

if k=l or normals < length then 
cio:=evalb(ttpl in Or); 
if cio=false then 

Or:=Or union transl(ttpl); Orl:=transl(ttp2);intcard:=nops(Or intersect Orl); 
if intcard > then NL[tnel+l]:=[n,mscofmt(ttpl),ttpl,sol[i]]; 
tnel:=tnel+l;else 

NL [tnel+ 1] : = [n,mscofmt (ttp 1 ) ,ttp 1 ,sol [i] ] ; 
NL [tnel+2] : = [n,mscofmt (ttp2) ,ttp2 ,sol [i] ] ; 
tnel:=tnel+2; 

fi; 
fi; 
fi; 

normals:=facesnew(ttp2);length:=checkrt(ttp2); 
if k=l or normals < length then 

ttp2:=tmttopq(mt[2]); 

cio:=evalb(ttp2 in Or); 
if cio=false then 

Or:=Or union transl(ttp2); Orl:=transl(ttpl);intcard:=nops(Or intersect Orl); 
if intcard > then NL[tnel+l]:=[n,mscofmt(ttpl),ttpl,sol[i]]; 
tnel:=tnel+l;else 

NL [tnel+ 1] : = [n,mscofmt (ttp 1 ) ,ttpl ,sol [i] ] ; 
NL[tnel+2]: = [n,mscofmt(ttp2),ttp2,sol[i]]; 
tnel:=tnel+2; 

fi; 
fi; 
fi; 

od; 
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[seq(L [i] ,i= 1 . .nel) ,seq(NL [j] ,j =nel+l . .tnel)] ; 
end: 

Once this list is computed for every n we can add up all the contributions for each tetrahedron 
in the list. A few procedures are used in the code of the main calculation called calculation. 

multbyfactorv:=proc(v,k) 
local w; 

w[l] :=v[l] *k;w[2] :=v[2] *k;w[3] :=v[3] *k; 
[w[l],w[2],w[3]]; 
end: 

multbyfactor.=pvoc(T,k) 
local i,NT,Q;NT:={}; 
Q : =convert (T , list ) ; 

for i from 1 to 4 do 

NT:=NT union multbyfactorv(Q[i],k); 

od;NT; 
end: 

addup:=proc(n,L) 

local i,j,m,mm,nt,k,T,Q,alpha,beta,gammma,d;nt:=0; 
m:=nops(L);k:=noor((n-l)/2); i:=l; 
while i<mdo 
if L[i][2] < n then 

mm:=floor(n/L[i][l]);T:=L[i][3];j:=l; 
while j < mm and L [i] [2] * j < n do 
Q : =multbyf actor ( T , j ) ; d : =L [i] [2] * j ; 

alpha:=nops(transl(Q));beta:=nops(intersalongEl(Q)); 
gammma:=nops(intersalongE2(Q));nt:=nt+f(n,d,alpha,beta,gammma); 

j:=j+i; 

od;fi; 

i:=i+l; 
od; 
nt/2; 
end: 

In the above procedure, addup, we get every irreducible from the list L together with all their 
appropriate dilations to compute their contribution in the cube [0, n] 3 using the formula given by 
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Theorem 13.21 The result is divided by two to match the sequence introduced in A103158. Finally, 
the list L is calculated in terms of n in two steps. First, construct the irreducible tetrahedrons 
using values of of mm = and nn = 1 in ExtendList and then we take care of the other possible 
values of mm and nn, making sure that we do not have duplicates by looking at the biggest d that 
shows up from each face. 

calculation:=pioc(n) local i,j,jl,j2,k,kl,L,x,xx,y,B,ii,T; 
x:=kvalues(n);y:=determined(n); xx:=floor((n+l)/2); 
i:=nops(x);L:=[];L:=ExtendList(l,L,0,l); 
for j from 1 to xx do 
L:=ExtendList(2*j+l,L,0,l); 
od; for j from 1 to i do 
k:=listofmn(x[j]);kl:=nops(k); 
for j2 from 1 to kl do 
for jl from 1 to nops(y[i]) do 
L:=ExtendList(y [j] [j 1] ,L,k[j2] [1] ,k[j2] [2] ) ; 
od; 
od; 

od; for ii from 1 to n do B[ii]:=addup(ii,L); print([ii,B[ii]]); od; 
T: = [seq(B[ii],ii=l..n)]; T; end: 

The result of the calculation(lOO) gives in less than a few hours of computations: [1, 1] [2, 9] [3, 
36] [4, 104] [5, 257] [6, 549] [7, 1058] [8, 1896] [9, 3199] [10, 5145] [11, 7926] [12, 11768] [13, 16967] 
[14, 23859] [15, 32846] [16, 44378] [17, 58977] [18, 77215] [19, 99684] [20, 126994] [21, 159963] [22, 
199443] [23, 246304] [24, 301702] [25, 366729] [26, 442587] [27, 530508] [28, 631820] [29, 748121] 
[30, 880941] [31, 1031930] [32, 1202984] [33, 1395927] [34, 1612655] [35, 1855676] [36, 2127122] 
[37, 2429577] [38, 2765531] [39, 3137480] [40, 3548434] [41, 4001071] [42, 4498685] [43, 5044606] 
[44, 5641892] [45, 6294195] [46, 7005191] [47, 7778912] [48, 8620242] [49, 9533105] [50, 10521999] 
[51, 11591474] [52, 12746562] [53, 13992107] [54, 15332971] [55, 16775590] [56, 18324372] [57, 
19985523] [58, 21765013] [59, 23668266] [60, 25702480] [61, 27873699] [62, 30188259] [63, 32655348] 
[64, 35281418] [65, 38074085] [66, 41040495] [67, 44188592] [68, 47525856] [69, 51061295] [70, 
54804647] [71, 58763604] [72, 62949850] [73, 67371219] [74, 72037311] [75, 76958126] [76, 82143618] 
[77, 87606245] [78, 93355379] [79, 99403446] [80, 105762770] [81, 112443331] [82, 119456581] [83, 
126814970] [84, 134532746] [85, 142621185] [86, 151093691] [87, 159964136] [88, 169245226] [89, 
178954039] [90, 189102295] [91, 199706864] [92, 210781424] [93, 222341631] [94, 234402515] [95, 
246978962] [96, 260093046] [97, 273757925] [98, 287989943] [99, 302809940] [100, 318235290] 
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Figure 6: The graph '"ffffi/j^ , 1 < n < 100. 
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d Finding the solutions a,b,c of a A 2+b A 2+c A 2=3d A 2 (0<=a<=b<=c). 

\()nly the solutions that satisfy gcd(a,b,c)=l are in the output. 

> abcsol : =proc (d) 
local i,j,k,m,u,x,y, sol , cd; sol :={ } ; 
for i from 1 to d do 

u:=[isolve(3*d A 2-i A 2=x A 2+y A 2) ] ;k:=nops(u) ; 
for j from 1 to k do 

if rhs (u[ j] [1] ) >=i and rhs (u [ j ] [2] ) >=i then 
cd:=gcd(gcd(i,rhs(u[j] [1])) ,rhs(u[j] [2])) ; 
if cd=l then sol : =sol union 
{sort([i,rhs(u[j] [1] ) , rhs (u [ j ] [2])])};fi; 
fi; 
od; 
od; 

convert (sol , list) ; 
end: 

-J Finding the two minimal regular tetrahedra given a,b,c such that 
a A 2+b A 2+c A 2=3d A 2 

The fourth point is the origin. We have included only two since all the other are 
going to be produced by cube transformations. 

A few other programs check if the output is what is suppose to be. 

> findpar :=proc (a,b,c) 
local i , j , sol ,mx ,nx, r , s ,my , ny , q,d,u , v, w, x, y , z ,mu , nu,mv, nv,ef , ns ,mz ,nz 
, mw , nw , oml , om2 , 1 , uu , t , Rl , R2 , f c ; q : =a A 2+b A 2 ; 

sol:=convert ( {isolve (2*q=x A 2+3*y A 2) } , list) ;ns :=nops (sol) ;d:=sqrt ( (a A 2+: 
A 2+c A 2)/3) ; 



and rau=f loor (mu) and nu=floor(nu) and mv=f loor (rav) and nv=f loor (nv) 
then 

u : =mu*m-nu*n ; v : =rav*m-nv*n ; w : =mw*m-nw*n ; 
x : =mx*m-nx*n ;y : =my*m-ny*n ; z : =mz*m-nz*n ; 
Rl:=[ (x+u-2*a)/3, (v+y-2*b)/3, (z+w-2*c)/3] ; 
R2:=[ (x+u+2*a)/3, (v+y+2*b)/3, (z+w+2*c)/3] ; 
fc:=subs (m=l,n=0 ,R1 [1] ) ; 
if f c=f loor (f c) then 

oml:=subs (m=l,n=0, [ [u,v,w] , [x,y,z] ,R1] ) ; else 
oml:=subs (m=l,n=0, [ [u,v,w] , [x,y,z] ,R2] ) ;fi; 
fc:=subs (m=l,n=l ,R1 [1] ) ; 
if fc=floor(fc) then 

om2 :=subs (m=l ,n=l , [ [u,v,w] , [x,y , z] ,R1] ) ; else 
om2 : =subs (m=l ,n=l , [ [u, v, w] , [x,y , z] ,R2] ) ;f i; 
ef:=l; fi;fi; 
od; 

oml , om2 ; 
end: 

> subtrv:=proc (U,V) 
local W; 

W[l] :=U[1]-V[1] ;W[2] :=U[2]-V[2] ;W[3] :=U[3]-V[3] ; 
[W[l] ,W[2] ,W[3]] ; 
end: 



jpj Translating the minimal tetrahedra to the positive quadrant. 

> tmttopq : =proc (T) 

local i,a,b,c,v,0,TR; 
a:=min(T[l] [1] ,T[2] [1] ,T[3] [1] ,0) ; 
b:=min(T[l] [2] ,T[2] [2] ,T[3] [2] ,0) ; 
c:=min(T[l] [3] ,T[2] [3] ,T[3] [3] ,0) ; 
O:=[0,0,0] ; 

v:=[a,b,c] ;TR:={subtrv(0,v) ,subtrv(T[l] ,v) ,subtrv(T[2] ,v) ,subtrv(T[3] , 

) }; 

end: 



.d Determining the number of tetrahedra in the orbit under the cube 
symmetries 



> orbitl : =proc (T) 

local i,k,Tl,a,b,c,x,T2,T3,T4,T5,T6,T7,T8,T9,T10,Tll,T12,T13,T14,T15,T16,T17,T18, 

T19,T20,T21,T22,T23,T24,S,Q,d,al,bl,cl; 

Q : =convert (T , list) ; 

d:=mscofmt(T) ; 

Tl :=T; 

T2:={seq([Q[k] [2] ,Q[k] [3] ,Q[k] [1] ] ,k=l . . 4) ) ; 
T3:={seq([Q[k] [1] ,Q[k] [3] ,Q[k] [2] ] ,k=l. .4) } ; 
T4:={seq([Q[k] [1] ,Q[k] [2] ,d-Q[k] [3]],k=l. .4) > ; 
T5:={seq([Q[k] [2] ,Q[k] [3] ,d-Q[k] [1] ] ,k=l. .4) > ; 
T6:={seq([Q[k] [1] ,Q[k] [3] ,d-Q[k] [2] ] ,k=l. .4) > ; 
T7:={seq([Q[k] [1] ,d-Q[k] [2] ,Q[k] [3] ] ,k=l. .4) >; 
T8:={seq([Q[k] [2] ,d-Q[k] [3] ,Q[k] [1] ] ,k=l. .4) >; 
T9:={seq([Q[k] [1] ,d-Q[k] [3] ,Q[k] [2] ] ,k=l. .4) > ; 
T10:={seq([d-Q[k] [1] ,Q[k] [2] ,Q[k] [3] ] ,k=1..4)}; 
Tll:={seq([d-Q[k] [2] ,Q[k] [3] ,Q[k] [1] ] ,k=l . . 4) } ; 
T12:={seq([d-Q[k] [1] ,Q[k] [3] ,Q[k] [2] ] ,k=l . . 4) } ; 
T13:={seq([Q[k] [1] ,d-Q[k] [2] ,d-Q[k] [3] ] ,k=l . . 4) } ; 
T14:={seq( [Q[k] [2] ,d-Q[k] [3] ,d-Q[k] [1] ] ,k=l. .4) } ; 
T15:={seq([Q[k] [1] ,d-Q[k] [3] ,d-Q[k] [2] ] ,k=l . . 4) } ; 
T16:={seq([d-Q[k] [l],d-Q[k] [2],Q[k] [3] ] ,k=l . . 4) } ; 
T17:={seq([d-Qtk] [2] ,d-Q[k] [3] ,Q[k] [1] ] ,k=l . .4) } ; 
T18:={seq([d-Q[k] [l],d-Q[k] [3],Q[k] [2] ] , k=l . . 4) } ; 
T19:={seq([d-Q[k] [1] ,Q[k] [2] ,d-Q[k] [3] ] ,k=l . . 4) } ; 
T20:={seq([d-Q[k] [2] ,Q[k] [3] ,d-Q[k] [1] ] , k=l . . 4) } ; 
T21:={seq([d-Q[k] [1] ,Q[k] [3] ,d-Q[k] [2] ] ,k=l. .4) } ; 
T22:={seq([d-Q[k] [1] ,d-Q[k] [2] ,d-Q[k] [3]],k=l. .4) }; 
T23:={seq([d-Q[k] [2],d-Q[k] [3],d-Q[k] [1] ] ,k=l . . 4) > ; 
T24:={seq([d-Q[k] [l],d-Q[k] [3],d-Q[k] [2] ] ,k=l . . 4) } ; 

S: = {T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11,T12,T13,T14,T15,T16,T17,T18,T19,T20,T21,T22, 
T23,T24} ; 

S; 

end: 

> orbit : =proc (T) 
local S,Q,T1; 

Q : =convert (T , list) ; 

Tl:={ [Q[l] [3] ,Q[1] [2] ,Q[1] [1] ] , [Q[2] [3] ,Q[2] [2] ,Q[2] [1] ] , [Q[3] [3] ,Q[3] [2] ,Q[3] [1] ] , 
[Q[4] [3] ,Q[4] [2] ,Q[4] [1] ] } ; 
S:=orbitl(T) union orbitl(Tl); 



addvect:=proc (T,v) 
local i,Q;Q:={}; 
for i from 1 to 4 do 
Q:=Q union { addvec (T [i] , v) } ; 
od; 

Q; 

end: 

transl : =proc (T) 

local S,Q,i,j,k,a2,b2,c2,a,b,c,d; 
Q:=convert (T, list) ; 

a:=max(Q[l] [1] ,Q[2] [1] ,Q[3] [1] ,Q[4] [1]) 
b:=max(Q[l] [2] ,Q[2] [2] ,Q[3] [2] ,Q[4] [2]) 
c:=max(Q[l] [3] ,Q[2] [3] ,Q[3] [3] ,Q[4] [3]) 
d: =max (a,b , c) ; 
a2 : =d-a ;b2 : =d-b ; c2 : =d-c ; 
S:=orbit(T) ; 
for i from to a2 do 
for j from to b2 do 
for k from to c2 do 

S:=S union orbit (addvect (T , [i , j , k] )) ; 
od; 

od; 
od; 
S; 

end: 

-J Intersection of the extended orbit with its translation along [0,0,1] 
(calculation of beta). 

> intersalongEl : =proc (T) 
local S,m,i,Sl,S2; 

S2 :=transl (T) ; S : =convert (S2 , list) ;m:=nops (S) ;Sl :={ } ; 
for i from 1 to m do 

S1:=S1 union {addvect (S [i] , [0 , , 1] )} ; 
od; 

S2 intersect SI; 



S1:=S1 union {addvect (S [i] , [0 , , 1] )} ; 
od; 

S3:={}; 

for i from 1 to m do 

S3:=S3 union {addvect (S [i] , [0 , 1 , 0] )} ; 
od; 

SI intersect S3 ; 
end: 

-J The function that gives the number of tetrahedra generated in an 
arbitrary cube. 

> f:=(n,d, alpha, beta,gammma) ->(n-d+l) A 3*alpha-3* (n-d) * (-d+l+n) A 2*beta+3* 
ammma* (n-d+1) * (n-d) "2 : 

> notetraincn:=proc (T,n) 
local d,x,y,z,w; 
d:=mscofmt (T) ; 

x:=nops (transl (T) ) ;y :=nops ( inter salongEl (T) ) ;w:=nops ( inter salongE2 (T) ) 
if n>=d then 

z :=f (n,d,x,y,w) ;else z:=0; 
fi; 
z ; 

end: 

List of minimal tetrahedra and the main programs. 

We need to keep for each minimal tetrahedra, the length of its sides (which is 
d\sqrt{2}), the dimension of the minimal cube containing it, the alpha, the beta, 
and the gamma numbers. 

> ExtendList:=proc(n,L) 
local 

i , sol , nsol , alpha , beta , gammma , nel , mt , ttpl , ttp2 , Orl , Or , tnel , NL , intcard , c 

o ; 

nel :=nops (L) ; 

sol:=abcsol (n) ; nsol :=nops (sol) ; 



NL[tnel+2] : = [n,mscofrat (ttp2) , ttp2 , sol [i] ] ; 
tnel : =tnel+2 ; 
fi; 
fi; 

cio : =evalb (ttp2 in Or); 
if cio=false then 

Or:=Or union transl (ttp2) ; Orl :=transl (ttpl) ;intcard:=nops 
intersect Orl) ; 

if intcard>0 then NL[tnel+l] :=[n,mscofmt(ttpl) ,ttpl,sol[i] ] 
tnel : =tnel+l ; else 

NL[tnel+l] : = [n,mscofmt (ttpl) ,ttpl,sol[i] ] ; 
NL[tnel+2] :=[n,mscofmt(ttp2) ,ttp2,sol[i] ] ; 
tnel : =tnel+2 ; 
fi; 
fi; 

od; 

[seq(L[i] , i=l . .nel) ,seq(NL[j] , j=nel+l. .tnel) ] ; 
end: 

> multbyfactorv: =proc (v,k) 
local w; 

w[l] :=v[l] *k;w[2] :=v[2] *k;w[3] :=v[3]*k; 
[w[l] ,w[2] ,w[3]] ; 
end: 

> multbyf actor : =proc (T , k) 
local i,NT,Q;NT:={ } ; 
Q:=convert (T, list) ; 
for i from 1 to 4 do 

NT : =NT union {multbyf actorv (Q [i] , k) } ; 

od ; NT ; 

end: 

> main:=proc (n,L) 

local i, j , m, mm, nt,k,T,Q, alpha, beta, gammma , d;nt : =0 ; 

m:=nops (L) ;k:=f loor ( (n-1) /2) ; 

#if 2*k+l=n then L:=ExtendList(n,L) : f i ; 

i : =1 ; 



od; 
nt/2 ; 
end: 

> calculation : =proc (n) 
local i,j,L; 
j :=f loor (n/2) ; 
L:=[] ;L:=ExtendList(l,L) ; 
for i from 1 to j+1 do 
L:=ExtendList(2*i+l,L) ; 
od; 

main (n, L) ; 
end: 

d Calculations of the list 

> calculation (1) ; calculation (2) ; calculation (3) ; calculation (4) ;calculatio 
(5) ; 

1 
9 

36 
104 
257 



> calculation (30) ; 

> calculation (31) ; 

> calculation (32) ; 

> calculation (33) ; 



SS0W1 



1031930 



12029&4 



1395927 



> calculation (34) ; 



